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Abstract — For the basic maximum likelihood estimating func- 
tion of the two parameters WeibuU distribution, a simple proof 
on its global monotonicity is given to ensure the existence and 
uniqueness of its solution. The boundary of the function's first- 
order derivative is defined based on its scale-free property. 
With a bounded derivative, the possible range of the root of 
this function can be determined. A novel root-flnding algorithm 
employing these established results is proposed accordingly, its 
convergence is proved analytically as well. Compared with other 
typical algorithms for this problem, the efficiency of the proposed 
algorithm is also demonstrated by numerical experiments. 

Index Terms — Two parameter WeibuU Distribution; Maximum 
Likelihood Estimation; Global Monotonicity; Scale-free Property. 



I. Introduction 

The Weibull distribution 1 1 1 is an important distribution in 
reliability and maintainability analysis?. The estimation on its 
parameters has been widely discussed, and there are several 
methodological categories for this parameter estimation issue 
10 lEI For example, the graphic methods ||5|, transcen- 
dental equation-solving method applying bifurcation algorithm 
l|41, maximum likelihood estimation (MLE) method [6] |7| 
lis). The graphic methods, such as Weibull probability plotting 
(WPP) |9|, are straight forward, but can not give a precise 
estimation. The transcendental equation-solving method has 
a closed form, it can avoid the computing expense in the 
iterative computation of the raw data, but in order to solve the 
transcendental equation, advanced mathematical techniques 
are required. 

In the MLE-based methods, as the basic estimating equation 
is not in closed form, it can be solved only numerically. 
There are several typical MLE-based methods for solving this 
equation, such as the secant method, the bisection method 
and the Newton-Raphson method. However, in both the secant 
method and the bisection method, the convergence rates are 
very low; in the Newton-Raphson method, it has to compute 
both the basic estimating function and its derivative [lOJ 
at each iterative step. In some cases, the Newton-Raphson 
method cannot ensure convergence fTT]. Furthermore, these 
above MLE-based methods require either initial values or trial 
computation of the estimated parameters. 
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In lfT2ll . the author claims that the existence and uniqueness 
of the solution of the basic estimating equation under MLE 
method cannot be assured. However, in |13| |14|, based on 
Cauchy-Schwarz inequality, the authors have given similar 
proofs of the existence and uniqueness on the solution of 
the MLE-based estimator In this paper, we present a straight 
forward proof by mathematics induction, which is much 
simpler than both the proofs given in ifTJl lfT4l . 

For the basic estimating function, we have proved the scale 
free property of its first order derivative, the boundary of 
this derivative can thus be defined. Moreover, with a bounded 
derivative, at each iterative step, the possible range of the root 
of the basic estimating function can be determined. We thus 
propose a novel MLE-based root-finding algorithm based on 
these properties, its computational efficiency and advantages 
are well demonstrated by numerical experiments. 

The remainder of this paper is organized as follows. In 
Section |II] a proof on the monotonicity of the basic estimating 
function is given. In Section |llll the scale-free property of the 
first order derivative of this function is proved, the boundaries 
of the function itself and its first order derivative are defined. 
The feasible range of its root thus can be determined. In 
Section IIVI by employing these proved results, a novel root- 
finding algorithm is hence designed. The convergence of the 
proposed algorithm is proved. In Section IIVI the performance 
of this proposed root-finding algorithm are demonstrated by 
numerical experiments. Conclusions are given in I VII 



II. The Global Monotonicity of the Basic 
Estimating Function 

A. The Basic Estimating Function for the Two Parameter 
Weibull Distribution 

The density function of the two parameter Weibull dis- 
tribution is: 

/(^) = (-)(^)'=e-(^)' (a; > 0, A: > 0, A > 0) (1) 

For a sampled data of n observations with the above EqnlT] 
as the applicable density function, its likelihood function is 



k 



Xj A 



(2) 
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By MLE method, the following equations can be obtained: 

9 In L n -J^ 1 ^ 

- + ^\nx,- TT 2^^^ Inx, = 0, 



dk 



91nL k 



1 " 



0. 



by eliminating A, we get 



Er=i 1^ _ i = i 



(3) 
(4) 

(5) 



by the above Eqn|5l we can get the value of k by related 
numerical algorithms. With k determined, A can be calculated 
by EqnUas 




A = = 



1 



(6) 



F(fc) = 



1 " 1 



(7) 



we can therefore calculate both k and A. Here k (A) refers to 
maximum likelihood estimators for parameter k (A). 

By Eqn|5] we can define the basic estimating function F{k) 

as 

Thus = here is defined as the basic estimating 

equation. By mathematics induction, with n > 1 and /c > 0, 
a simple proof on the global monotonicity of F(fc), i.e., 
> can be given in the following section. 

B. Proof on the global monotonicity of F{k) 
Since 

dF{k) _ 1 

fc2 + 



Fig. 1. The boundaries of F(k): < F{k) < + Ci. Here Ci 



III. Boundaries of the Basic Estimating Function 
AND its First-order derivative 

A. Boundaries of the Basic Estimating Function F{k) 
Let 



G{k) 



as 



G(0) = 



■El- 



(13) 



G(+oo) = iE:ri'ln^>0 
also because 

dG{k) 



dk 



dk 



Pi{x,,n,k){Y,x'iy^ >Q, (14) 



i—1 i—1 i—1 i—1 



(8) 



dFjk) 



to prove " g^"^ > 0, we need only to prove that for fc > and 
any Xi G i?"*", 

n n n 

Pi{xi,n,k) = ^xfln^x^^x*^ - (^xf Ina;^)^ > 



thus Vfc > 



here 



n— 1 

1 , a;?: 



i=l 



Ci = G(+oo) = - V In 

77 Z ^ 



for n = 1 

Pi{x^, l,k) = x'^ln^xi ■ x'l - xl'' ■ In^xi = 
for n — 2 

Pi{x„2,k) ^x'lx^{\iiX2 -Inxi)^ > 
for ri = TO — 1 (m > 3, TO G N^), suppose 
Pi{x„m- l,k)>0 

then for n = m 

Pi{xi,m, k) — Pi{xi,m — l,k) + 



n ^ Xi 



(15) 



(16) 



(9) 
(10) 
(11) 



1=1 

Therefor ^^g^ > 0, and function F{k) is global monotonic. 
With the global monotonicity of F{k), the existence and 
uniqueness of the root of F(k) = can be assured. 



we define the lower boundary of F{k) is curve: F2{k) = 
— the upper boundary of F{k) is curve: Fi{k) — — -i + Ci. 
Then the boundary curves of F{k) can be seen in FiglT] 



B. Boundaries of the First Order derivative: ^^q^^ 
Let 

^ n n 

H{n, x„ fc) = ] XI ^^'lii^^;^ X x,^ 

'^i=l i=l 

n n N n 

-J^xflnx.^x.^lnx, (Xxf)-2 (17) 

i=l i=l i=l 

Theorem A: Vn e A^+,Va;, > 0, and Vfc > 0, G 



[fc-2,fc-2 + ln2(^-^- 
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PROOF: For a certain Xj > satisfying ln(Xj ■ Xi) > 0, 

H(n, Xi, k) — Hin, Xj ■ Xi,k) 

< max{77(n, Xj ■ Xi, k)} 



In- 



ln(A^ ■ x^^^ ■ x^^J 



(18) 



Since to satisfy 



ln(Aj ■Xi)>0 

ln(A2 •a;,^ ■ x^^J > 
the minimum value of Xj is x^^^^, i.e., Aj G [a;^^^^ 

let A = 7 • (7 > 1), 



max {H{n, X ■ Xi,k)} = ) 111(7 — 

fc>0,7>l 

therefore the supermum value of H{n, Xi, k) is 
sup{H{n,Xi,k)} — min{ max {iJ(?T,, A • x^, A:)}} 

a;i>0,neA''"'', 
fe>0,7>l 

= mini ln(^^)[ln72 + ln(^^)l 
= ln^(^) 



-) (19) 



Hence 



H{n, Xi, k) < ln^( ""^ ) 



(20) 



It is proved in Section lTl-BI the infimum value of H{n, Xxi, k) 
is: 

M{H{n,x„k)} ^0 (21) 



with 



OF(fc) _ 
dk ~ 



H{n, Xi, k) 



therefore 



1 < ^ <C, + - 
k^ dk k"^ 



(22) 



here C2 -ln2(f^). □ 

With the boundary 
range of F{k) = can be determined 



With the boundary of defined, the feasible solution 



C. Possible Range of the Final Solution of Equation F{k) — 
The basic estimating function F{k) can be written as 

F{k) = I F'{T)dT 
Jko 

According to Theorem A, 

IL ^-'dr < F' ir)dr < [C, + r-^]dr (fc > ko) 

IL r-'dr > J' F'{r)dr > j'[C,+r-^]dr (k < ^o) 




FL(k,ko) 



Fig. 2. The two boundary curves of F{k): Fi^ik, ko) and Fu{k, ko), which 
can determine the feasible range of the solution of F{k) = 0. 



thus 

FL{k, ko) < F{k) < Fu{k, ko) (fc > fco) 

FL{k, ko) > F{k) > Fu{k, ko) (fc < fco) 

Here Fu{k,ko) and ^^(fc, fco) can be seen the boundary 
curves of F{k) at point (fco, F(fco)) (see Fig|2]l. With knowing 
F(fco), the two boundary curves of F{k) can be determined 

as, 

Fi(fc,fco) = fco-'+F(fco)- i 

Fu{k,ko) = Csfc - i + (F(fco) - Cafco + ^) 

We can see that both Fjjik, ko) and FL{k, ko) are algebraic 
functions with simple mathematical forms. 

Suppose Fjj{k,ko) cuts fc-axe at point fcj^, FL{k,ko) cuts 
fc-axe at point fc^, in Fig|2] it is obvious that the root of 
equation F(k, ko) = must exist in [fc^, fc^*,] (if fc° < k^) or 
in [fcO,fcO] (if fc^ >fcO). 

A detailed algorithm employing the properties in Section 
nil- Al IIII-B I and llll-CI will be proposed in the following section, 
and a proof on its convergence will be demonstrated as well. 



IV. Root-finding Algorithm Design 

To numerically get the solution of F{k) = 0, an intuitive 
idea is to calculate F{kl) at fci — ^(fc° + fc|j), to further 
find out a narrower interval [fc^^, k^] or [fcj^, kjj]. By the same 
procedure, at fci+i = ^{k\ + klj){i — > +oo), the final solution 
of F{k) = can be iteratively approximated. 

From the mathematical form of the basic estimating function 
F{k), we can see that for large n, k and Xi, the computational 
complexity for F{k) is very high, thus the calculation time 
of F{k) is an important index in measuring the efficiency of 
a root-finding algorithm for F{k) = 0. A novel root-finding 
algorithm has been designed based on these properties proved 
in Section |ll] , |III] and UlLCl 



let denote 



^L(fc,fco)=/, 



k 2 



dr 



Fu{k,ko)= Sl[C2+T-^]di 



A. Convergence of the Bounded derivative Algorithm 

As stated before, a straight forward method employing this 
idea is to calculate the final solution iteratively. The feasibility 
of this method can be ensured by the following theorems. 
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Casell: F(ki+i)<0 



' Fl(K, Ki) 



(ki,F(ki)) 



Fig. 3. Estimation on the convergence rate of tlie proposed algorithm with 
knowing F{ki) and i^(fci+i), here Fi^{k,ki) and Fu{k,ki) refer to the 
boundary functions of F{k) at point {ki, F{ki)), ki^i = ^(k]^ + k^). 



Lemma. For any point {ki,F{ki)) of function F{k), the 
existence and uniqueness of the roots of Fij{k,ki) — 
and F]^{k,ki) = can be assured, therefore the interval 
[k\^,k\j\ (i > 1) that covers the root of F{k) — always 
exists. 

PROOF: (1) By Inequality [H we know that VA: > 0, 

F{k) > -i. The root of equation FL(fc,A;o) = + F{ko)- 



k°. 



1 



F{ko) - i-k^') 



> 



(23) 



which is reasonable. 

(2) Since Fij{k, kg) = can be rewritten as 

1 



C2k' + iF{ko) - Cafco + — )fc -1 = (24) 

fco 



Since here 



A = B^- 4AC = {F{ko) - C2ko + + AC2 > 
C = -1< 



it is obvious that Fij{k,ko) = always has two roots with 
different signatures. To satisfy fc > 0, only the positive root 
should be preserved. 

Therefore, the existence of the interval [k]^, k'jj] {i > I) can 
be assured. □ 

Theorem B. Let [k]^ , A:^] denote the feasible interval of 
the solution of F{k) = at point {ki,F[ki)), 
denote the feasible interval at point (fc^+i, ^(fci+i)). 



with fc 



+1 2 



k\j), then the convergence rate 

i. Here fc)+i = max(fci^, fc)+i). 



k\J' = min(fc^, fc[+^), i e iV^ 



PROOF: With knowing point {ki,F{h)) {i > 1), we can 
deduce the boundary curves: FL{k,ki) and Fij{k^ki), and 
the related feasible interval [k^^^ , k^^^] can be determined 
accordingly (as shown in Fig|3]l. 

Since if F(fci+i) ~ 0, then K is the root we are looking 
for. Thus here there are only two situations for F{ki+i): 
F(fc,+i) > and F(h+i) < 0. 

Case I: F{h+i) > 



Fig. 4. The flow chart for the combined algorithm. Here Si is a scale 
threshold, 62 is the approximating precision or the computing halt criterion. 



With knowing F{ki+i) [i > 1), the boundary curves 
FL{k, ki+i) and Fu{k, fcj+i) at point {ki+i, F{ki+i)) can be 
determined, and so is the related feasible interval [fc)^^^, fcj/^]- 

As in this case, by Figl3] it is obvious that k^^ < Ki^i, 
then the convergence rate 



7 



< 



(25) 



Case IL F{h+i) < 

Similarly, as shown in Figl3] in this case, it is obvious that 
k^]^^ > Ki+i, then the convergence rate 



7 ^ 



< 



< 



fc. 



kl 



Thus the proposed algorithm is convergence, and its con- 



vergence rate 7 < 



□ 



B. Improved Algorithm Combining the Secant Method and the 
Bounded derivative Method 

It is obvious that at large scale the bounded derivative 
algorithm can converge rapidly, at least at the same rate as 
the bisection method. While in a small scale, especially in 
the linearizable neighborhood around the final solution point 
of the basic estimating equation, the secant method method 
has a better convergence rate. Therefore it is intuitive to 
combine both the two methods together. The flow chart for 
the combined algorithm can be seen in Fig ID 

As shown in the flow chart in FiglH in the combined 
algorithm, at a large scale, at each iterative step, we calculate 
F(ki) one time, then with only basic algebraic calculation, we 
find the possible range of the solution, and then narrow down 
it to a smaller range at next step. When the possible range 
of the root comes to a small scale, where the linearization 
part of the basic estimating function dominates, the combined 
algorithm will switch to the secant method, since it has a better 
performance in this case. 
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Moreover, compared with the secant method, the bisection 
method and the Newton-Raphson method, which requires 
initial values of F{k), no initial value is required in the 
combined method. 

V. Numerical Examples 

A. General Cases 

With MATLAB, 1000 times of numerical experiments have 
been taken to study the performance of this proposed algo- 
rithm. In these simulations, the shape parameter of Weibull 
distribution k ^ [/(0,40], the scale parameter A ~ f/(0,40], 
the number of data N ^ C/[2, 1000], here U refers to Uniform 
Distribution. 

Compared with the secant method, the bisection method 
as well as the bi-secant method (which has combined the 
secant and the bisection method), the computing complexity 
(here refers to the calculation time of F{k)) of each method 
under various approximating precision e can be seen in the 
following Table I. 

Table I: Averaged calculation times of F{k) 



Appro - 


Secant 


Bisection 


Bi-secant 


Proposed 


precision e 


method 


method|6| 


Method 


method 


10-^ 


9.42 


5.21 


3.12 


1.68 


10-^ 


30.22 


6.33 


4.29 


2.58 


10'^ 


88.5 


11.22 


9.39 


4.28 


10-* 


130.6 


14.19 


11.51 


5.05 



B. A Constructed Case 

The advantage of the proposed algorithm can also be 
verified by a concrete case, which is generated from a Weibull 
Distribution. The sampled data is in Table II, the comparison 
of the performance of different methods can be seen in Table 
III. 

Table ILSampled Data from a Weibull Distribution 



2.6144 


4.1834 


4.3258 


4.3496 


4.3740 


4.4006 


3.2073 


4.2573 


4.3273 


4.3544 


4.3828 


4.4051 


3.9800 


4.2884 


4.3334 


4.3646 


4.3873 


4.4123 


4.1767 


4.3150 


4.3403 


4.3698 


4.3959 


4.4194 


4.4317 


4.4919 


4.4448 


4.5082 


4.4623 


4.5439 


4.4756 


4.5715 











Table III: Calculation times of F{k) of various methods 



scale-free property of this function's first order derivative, the 
possible range of its solution at each iterative step can be 
determined. Based on these properties, a novel root-finding 
algorithm is proposed in this paper. Its efficiency has been 
demonstrated by numerical experiments. 

For future work, this method can be extended to type I and 
type II data censoring cases. 
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Appro- 
precision € 


Secant 
method 


Bisection 
method 


Bi-Secant 
Method 


Proposed 
method 




4 


4 


3 


1 


IQ-^ 


79 


6 


3 


1 


IQ-^ 


178 


10 


3 


2 


10-4 


271 


13 


21 


3 


10-" 


459 


20 


22 


4 


10-^" 


833 


33 


23 


5 


10-^4 


1206 


46 


24 
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VI. Conclusion 

In this paper, to assure the uniqueness and existence of the 
solution of the basic estimating function, we have proved its 
global monotonicity in a very simple way. With proving the 



